Naszym zadaniem jest przeanalizowanie zbioru danych medycznych oraz stworzenia modelu klasyfikującego potencjalnego pacjenta ze względu na fakt wystąpienia udaru mózgu.
Udar mózgu to zespół objawów klinicznych związanych z nagłym wystąpieniem ogniskowego lub uogólnionego zaburzenia czynności mózgu, powstały w wyniku zaburzenia krążenia mózgowego i utrzymującego się ponad 24 godziny.
import pandas as pd
import numpy as np
import joblib
from plotly.subplots import make_subplots
import plotly.graph_objects as go
import plotly.express as px
from matplotlib import pyplot as plt
import seaborn as sns
from sklearn import preprocessing
from sklearn.utils import shuffle
from sklearn.model_selection import train_test_split
from imblearn.over_sampling import SMOTE
from imblearn.over_sampling import SMOTENC
from imblearn.over_sampling import ADASYN
from imblearn.combine import SMOTETomek
from imblearn.combine import SMOTEENN
from imblearn.under_sampling import EditedNearestNeighbours
from sklearn.pipeline import Pipeline
from sklearn import metrics
from sklearn.model_selection import GridSearchCV
from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
from sklearn.feature_selection import SequentialFeatureSelector
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import RFE
from sklearn.model_selection import cross_validate
from sklearn.feature_selection import SelectFromModel
from yellowbrick import classifier
from yellowbrick.classifier import ClassificationReport
from yellowbrick.classifier.rocauc import roc_auc
import warnings
warnings.filterwarnings("ignore")
def plotting(data: pd.DataFrame, x: str, type: str, y=None, cl = 'stroke'):
"""
Function for plotting data in specific way
Args:
data - dataset to plot
x - column for x-axis
y - column for y-axis
type - type of plot
values: {'hist', 'box', 'hist_box', 'class_hist', 'class_box', 'class_hist_box', 'scatter', 'heatmap', 'lineplot'}
cl - class
Returns:
Plot in plotly
"""
def plot_class_hist():
"""
Plots histogram with stroke classes in plotly
"""
fig = go.Figure()
fig.add_trace(go.Histogram(x = data.loc[data[cl]==0, x], name='No '+cl))
fig.add_trace(go.Histogram(x = data.loc[data[cl]==1, x], name=cl))
fig.update_layout(barmode='overlay')
fig.update_traces(opacity=0.75)
return fig
def plot_class_box():
"""
Plots boxplot with stroke classes in plotly
"""
fig = go.Figure()
fig.add_trace(go.Box(y = data.loc[data[cl]==0, x], name='No '+cl))
fig.add_trace(go.Box(y = data.loc[data[cl]==1, x], name=cl))
fig.update_layout(barmode='overlay')
fig.update_traces(opacity=0.75)
return fig
if type=='hist':
fig = px.histogram(data, x=x)
fig.show()
elif type=='box':
fig = px.box(data, x=x)
fig.show()
elif type=='scatter':
data['temp'] = data[cl].astype('str')
fig = px.scatter(data, x=x, y=y, color='temp')
fig.update_layout(title_text="Scatterplot for " + x)
fig.update_traces(marker_size=10)
fig.show()
elif type=='heatmap':
corr = data.corr()
f, ax = plt.subplots(figsize=(16, 12))
mask = np.triu(np.ones_like(corr, dtype=bool))
cmap = sns.diverging_palette(230, 20, as_cmap=True)
sns.heatmap(corr, annot=True, mask = mask, cmap=cmap)
elif type=='hist_box':
fig = make_subplots(rows=1, cols=2)
fig.add_trace(go.Histogram(x = data[x], name=x), row=1,col=1)
fig.add_trace(go.Box(y = data[x], name=x), row=1, col=2)
fig.update_layout(barmode="overlay")
fig.update_traces(opacity=0.7, row=1, col=1)
fig.update_layout(title_text="Histogram and Boxplot for " + x)
fig.show()
elif type=='class_hist':
fig = plot_class_hist()
fig.show()
elif type=='class_box':
fig = plot_class_box()
fig.show()
elif type=='class_hist_box':
fig = make_subplots(rows=1, cols=2)
fig.add_trace(go.Histogram(x = data.loc[data[cl]==0, x], name='No '+cl), row=1,col=1)
fig.add_trace(go.Histogram(x = data.loc[data[cl]==1, x], name=cl), row=1, col=1)
fig.add_trace(go.Box(y = data.loc[data[cl]==0, x], name='No '+cl), row=1, col=2)
fig.add_trace(go.Box(y = data.loc[data[cl]==1, x], name=cl), row=1, col=2)
fig.update_layout(barmode="overlay")
fig.update_traces(opacity=0.7, row=1, col=1)
fig.update_layout(title_text="Histogram and Boxplot with classes for " + x)
fig.show()
elif type=='lineplot':
fig = px.line(data, x=x, y=y, color=cl)
fig.show()
else:
raise NotImplementedError('Not implemented!')
return 0
def iqr_outliers(df: pd.DataFrame, col: str):
"""
Function for managin outliers in given column from dataframe using IQR method.
In descriptive statistics, the interquartile range (IQR) is a measure of statistical dispersion, which is the spread of the data.
It is defined as the difference between the 75th and 25th percentiles of the data.
Args:
df - dataset to remove outliers
col - column for removing outliers
Returns:
Whole dataframe with specified column cleaned from outliers using IQR method
"""
Q1 = df[col].quantile(0.25)
Q3 = df[col].quantile(0.75)
IQR = Q3 - Q1
df.loc[df[col] <= Q1 - 1.5*IQR, col] = df[col].quantile(0.1)
df.loc[df[col] >= Q3 + 1.5*IQR, col] = df[col].quantile(0.9)
return df
def balance_data(X, y, algorithm, cat_features, strategy=1, seed = 123, k_neighbors = 5, n_neighbors=3):
"""
Function for balancing dataset using specified algorithm and sampling strategy.
Args:
X - matrix of features
y - vector of target variable observations
algorithm - algorithm for balancing dataset
values: {'SMOTE', 'SMOTETOMEK', 'SMOTEENN'}
strategy - sampling information to resample the data set
Optional **kwargs:
k_neighbors - number of neighbours for SMOTE algorithm
n_neighbors - number of neighbours for ADASYN and ENN algorithms
Returns:
Whole balanced dataframe
"""
ft = []
for feat in cat_features:
ft.append(X.columns.get_loc(feat))
#sm = SMOTE(sampling_strategy = strategy, k_neighbors = k_neighbors)
sm = SMOTENC(random_state=seed, categorical_features=ft, sampling_strategy = strategy, k_neighbors=k_neighbors)
if algorithm == 'SMOTE':
bs = sm
elif algorithm == 'SMOTETOMEK':
bs = SMOTETomek(random_state = seed, sampling_strategy = strategy, smote= sm)
elif algorithm == 'SMOTEENN':
enn = EditedNearestNeighbours(n_neighbors = n_neighbors)
bs = SMOTEENN(sampling_strategy = strategy, smote=sm, enn = enn)
else:
raise NotImplementedError('Not implemented!')
X_bs, y_bs = bs.fit_resample(X, y)
df_bs = X_bs.copy()
df_bs['stroke'] = y_bs
X_bs = X_bs.drop(columns='index', axis=1)
print("Liczba wszystkich obserwacji: ", len(y_bs))
print("\nLiczba obserwacji z klasy mniejszościowej: ", len(y_bs[y_bs==1]))
print("\nProcent klasy mniejszościowej do całości po zbalansowaniu: ", len(y_bs[y_bs == 1])/len(y_bs))
return df_bs, X_bs, y_bs
#Stwórz klase z doborem parametru thresholda p
class LogisticRegressionWithThreshold(LogisticRegression):
def predict(self, X, threshold=None):
if threshold == None: # If no threshold passed in, simply call the base class predict, effectively threshold=0.5
return LogisticRegression.predict(self, X)
else:
y_scores = LogisticRegression.predict_proba(self, X)[:, 1]
y_pred_with_threshold = (y_scores >= threshold).astype(int)
return y_pred_with_threshold
class SupportVectorMachineWithThreshold(SVC):
def predict(self, X, threshold=None):
if threshold == None: # If no threshold passed in, simply call the base class predict, effectively threshold=0.5
return SVC.predict(self, X)
else:
y_scores = SVC.predict_proba(self, X)[:, 1]
y_pred_with_threshold = (y_scores >= threshold).astype(int)
return y_pred_with_threshold
class RandomForestClassifierWithThreshold(RandomForestClassifier):
def predict(self, X, threshold=None):
if threshold == None: # If no threshold passed in, simply call the base class predict, effectively threshold=0.5
return RandomForestClassifier.predict(self, X)
else:
y_scores = RandomForestClassifier.predict_proba(self, X)[:, 1]
y_pred_with_threshold = (y_scores >= threshold).astype(int)
return y_pred_with_threshold
def preprocess_data(X_train, X_test, columns_to_standardize, cat_features, scaler=True, ohe=True, test_only=False, sc=None):
"""
Function preprocesses data by standardizing numerical features and one hot encoding categorical features.
Args:
X_train - train dataset of features
X_test - test dataset of features
columns_to_standardize - columns to standardize
cat_features - column for OHE
scaler - boolean value specifies if standardizing should be performed
ohe - boolean value specifies if OHE should be performed
test_only - flag if only test set should be preprocessed
sc - scaler
Returns
X
Preprocessed train data
X_t
Preprocessed test data
"""
def ohe(df: pd.DataFrame, cat_features):
df_ohe = df[cat_features].astype('str')
dummies = pd.get_dummies(df_ohe)
df = df.drop(columns=cat_features, axis=1)
df = pd.merge(df, dummies, left_index=True, right_index=True)
return df
if test_only==True:
X_t = X_test.copy()
if scaler == True:
X_t[columns_to_standardize] = sc.transform(X_t[columns_to_standardize])
if ohe == True:
X_t = ohe(X_t, cat_features)
return X_t
else:
X = X_train.copy()
X_t = X_test.copy()
sc = preprocessing.StandardScaler()
if scaler == True:
X[columns_to_standardize] = sc.fit_transform(X_train[columns_to_standardize])
X_t[columns_to_standardize] = sc.transform(X_t[columns_to_standardize])
if ohe == True:
X = ohe(X, cat_features)
X_t = ohe(X_t, cat_features)
return X, X_t, sc
def fit_pred_score(X_train, y_train, X_test, y_test, model, model_name, dataset_name, scaler=True, ohe = False, columns_to_standardize = None, cat_features = None, visualize_train = True, visualize_test = True, show_coefficients = False, threshold = None):
"""
Function for training and testing given classifier.
It plots classification reports for train and test data and might display coefficients of predictors.
Returns report in dataframe format for given model and dataset with metrics from classification report.
Args:
X_train, y_train, X_test, y_test - training and testing data with associated labels
model - model to train and test
model_name - name of the model for the report
dataset_name - dataset name used for the report
scaler - if features should be scaled
columns_to_standardize - columns which should be scaled if scaler = True
visualize_train - if classification report for training data should be visualized
visualize_test - if classification report for testing data should be visualized
show_coefficients - if coefficients of features should be displayed
threshold - value of 'p' for Logistic Regression model
Returns:
report - dataframe with metrics from classification report with associated dataset and model name
"""
X, X_t, standard_scaler = preprocess_data(X_train, X_test, columns_to_standardize = columns_to_standardize, cat_features = cat_features, scaler = scaler, ohe = ohe)
model.fit(X, y_train)
if threshold == None:
y_pred_train = model.predict(X)
y_pred = model.predict(X_t)
else:
y_pred_train = model.predict(X, threshold = threshold)
y_pred = model.predict(X_t, threshold = threshold)
classes=['No stroke', 'Stroke']
#plot classification report for training data
if visualize_train == True:
visualizer = ClassificationReport(
model, classes=classes)
visualizer.fit(X, y_train)
visualizer.score(X, y_train)
visualizer.show()
if visualize_test == True:
#plot classification report for test data
classifier.classification_report(model, X, y_train, X_t, y_test, classes=classes)
if show_coefficients == True:
coef_dict = {}
for coef, feat in zip(model.coef_[0,:],X.columns):
print(feat, " = ", coef)
report = pd.DataFrame()
for t, y, y_approx in zip(['train','test'], [y_train, y_test], [y_pred_train, y_pred]):
results = pd.DataFrame(metrics.classification_report(y, y_approx, target_names = classes, output_dict=True))
results['metric'] = results.index
results.reset_index(inplace=True, drop=True)
results['dataset_name'] = dataset_name
results['model_name'] = model_name
results['split'] = t
report = pd.concat([report, results])
return report
def transform_report_to_plot(dt : pd.DataFrame, x : str, cl = 'metric', train_test = False):
"""
Function transforms report dataset for plotting in a way that for F1-score macro avg value is taken, for precision the values for no stroke and for recall the values for stroke.
Args:
dt - filtered dataset
x - axis
cl - class to display
train_test - flag which specifies if data should be transformed to differences between train and test
Returns:
df - transformed report dataframe
"""
if train_test == True:
prefix = 'train_test_diff'
else:
prefix = ''
#cols_to_select
f1_s = dt.loc[(dt.metric == 'f1-score'),[prefix+'macro avg', cl, x]]
precision_s = dt.loc[(report.metric == 'precision'),[prefix+'Stroke', cl, x]]
precision_s[cl] = "Precision (Stroke)"
recall_s = dt.loc[(dt.metric == 'recall'),[prefix+'Stroke', cl, x]]
recall_s[cl] = "Recall (Stroke)"
precision_ns = dt.loc[(report.metric == 'precision'),[prefix+'No stroke', cl, x]]
precision_ns[cl] = "Precision (No stroke)"
recall_ns = dt.loc[(dt.metric == 'recall'),[prefix+'No stroke', cl, x]]
recall_ns[cl] = "Recall (No stroke)"
columns = ['value', cl, '']
for d in [f1_s, precision_s, recall_s, precision_ns, recall_ns]:
d.columns = columns
df = pd.concat([f1_s, precision_s, recall_s, precision_ns, recall_ns])
return df
def prepare_diff_dataset(df : pd.DataFrame(), by: str, cols_for_diff = ['No stroke', 'Stroke', 'accuracy', 'macro avg', 'weighted avg']):
"""
Takes the report dataset and returns dataset with calculated differences in metric between test and train by specified grouping.
Args:
df - dataset
by - specified grouping
cols_for_diff - cols to calculate differences
Returns:
df_diff
Dataset with calculated differences
"""
df_diff = df.copy()
for col in cols_for_diff:
df_diff['train_test_diff'+col] = None
for ds in df_diff[by].unique():
for m in df_diff.metric.unique():
filter = (df_diff[by] == ds) & (df_diff.metric == m)
df_diff.loc[filter, 'train_test_diff'+col] = np.abs(df_diff.loc[filter, col].diff())
return df_diff
url = 'https://github.com/maju116/Statistics_II_GUT/raw/main/PROJECT/healthcare-dataset-stroke-data.csv'
df = pd.read_csv(url)
df.head(5)
| id | gender | age | hypertension | heart_disease | ever_married | work_type | Residence_type | avg_glucose_level | bmi | smoking_status | stroke | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 9046 | Male | 67.0 | 0 | 1 | Yes | Private | Urban | 228.69 | 36.6 | formerly smoked | 1 |
| 1 | 51676 | Female | 61.0 | 0 | 0 | Yes | Self-employed | Rural | 202.21 | NaN | never smoked | 1 |
| 2 | 31112 | Male | 80.0 | 0 | 1 | Yes | Private | Rural | 105.92 | 32.5 | never smoked | 1 |
| 3 | 60182 | Female | 49.0 | 0 | 0 | Yes | Private | Urban | 171.23 | 34.4 | smokes | 1 |
| 4 | 1665 | Female | 79.0 | 1 | 0 | Yes | Self-employed | Rural | 174.12 | 24.0 | never smoked | 1 |
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 5110 entries, 0 to 5109 Data columns (total 12 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 id 5110 non-null int64 1 gender 5110 non-null object 2 age 5110 non-null float64 3 hypertension 5110 non-null int64 4 heart_disease 5110 non-null int64 5 ever_married 5110 non-null object 6 work_type 5110 non-null object 7 Residence_type 5110 non-null object 8 avg_glucose_level 5110 non-null float64 9 bmi 4909 non-null float64 10 smoking_status 5110 non-null object 11 stroke 5110 non-null int64 dtypes: float64(3), int64(4), object(5) memory usage: 479.2+ KB
df.describe().transpose()
| count | mean | std | min | 25% | 50% | 75% | max | |
|---|---|---|---|---|---|---|---|---|
| id | 5110.0 | 36517.829354 | 21161.721625 | 67.00 | 17741.250 | 36932.000 | 54682.00 | 72940.00 |
| age | 5110.0 | 43.226614 | 22.612647 | 0.08 | 25.000 | 45.000 | 61.00 | 82.00 |
| hypertension | 5110.0 | 0.097456 | 0.296607 | 0.00 | 0.000 | 0.000 | 0.00 | 1.00 |
| heart_disease | 5110.0 | 0.054012 | 0.226063 | 0.00 | 0.000 | 0.000 | 0.00 | 1.00 |
| avg_glucose_level | 5110.0 | 106.147677 | 45.283560 | 55.12 | 77.245 | 91.885 | 114.09 | 271.74 |
| bmi | 4909.0 | 28.893237 | 7.854067 | 10.30 | 23.500 | 28.100 | 33.10 | 97.60 |
| stroke | 5110.0 | 0.048728 | 0.215320 | 0.00 | 0.000 | 0.000 | 0.00 | 1.00 |
Mamy łącznie 12 kolumn po 5110 obserwacji z czego jedna z nich stanowi ('stroke') zmienną celu i będziemy starać się ja zamodelować. Jest to zmienna binarna informująca czy dany pacjent miał udar mózgu czy też nie zatem jest to typowe zadanie klasyfikujące. Oprócz tego jest pozostałych 11 zmiennych:
Widzimy, że mamy 211 braków danych dla zmiennej bmi, którymi zajmiemy się w dalszej części projektu.
df.stroke.value_counts()
0 4861 1 249 Name: stroke, dtype: int64
plotting(df, 'stroke', type = 'hist')
0
Widzimy, że nasza zmienna celu jest niezbalansowana gdzie klasą większościową są pacjenci bez udaru mózgu. Musimy mieć to na uwadze i w przypadku pozostawienia takiego stanu rzeczy nie możemy skorzystać z metryki accuracy, a powinniśmy skorzystać z innej - przykładowo średniej harmonicznej między precision a recall czyli f1-score. Niezbalansowany zbiór danych jest jednak dość problemtyczny, gdyż klasyfikator ma tendencję skupiania się na predykcji klasy większościowej.
Istnieją sposoby na zbalansowanie zbioru danych, które są częścią pre-processingu czyli wstępnego przetwarzania danych przed modelowaniem. Do jednych z bardziej popularnych metod należą:
W przypadku undersamplingu możemy pozbyć się istotnych informacyjnie obserwacji. Dodatkowo nasz zbiór danych jest niewielki, więc nie chcielibyśmy się ograniczyć do zaledwie 249 obserwacji. Natomiast oversampling może doprowadzić do nadmiernego dopasowania modelu do naszych danych. Istnieje też podejście mieszane wykorzystujące zarówno undersampling, jak też oversampling ale ciężko jest ustalić optymalną proporcję między podejściami. Spróbujemy sprawdzić empirycznie najlepszą metodę i się do niej ograniczyć oczywiście opierając się nie o zwykłe próbkowanie, a dostępne i znane algorytmu balansowania zbiorów danych.
print("Liczba zduplikowanych wartości ID pacjentów: ", len(df[df.id.duplicated()]))
Liczba zduplikowanych wartości ID pacjentów: 0
W naszym zbiorze danych każdy pacjent jest wpisany unikalnie, tak więc zmienna ID raczej do niczego więcej nam się nie przyda i będziemy mogli ją pozostawić w dalszych etapach.
plotting(data = df, x='gender', type='class_hist')
0
df.gender.value_counts()
Female 2994 Male 2115 Other 1 Name: gender, dtype: int64
W zmiennej gender występują 3 wartości (male, female, other). Wartość 'Other' występuje tylko jeden raz, więc można go potraktować jako brak danych. W zależności od wartości pozostałych zmiennej dla tej obserwacji albo zastąpimy ją wartością male lub female (jako tą która najczęściej występuje) albo usuniemy.
plotting(data = df, x='age', type='hist_box')
0
plotting(data = df, x='age', type='class_hist_box')
0
Możemy zaobserwować, że udar mózgu w naszym zbiorze mają głównie osoby w podeszłym wieku (60-80 lat). Natomiast osoby bez udaru rozkładają się dość równomiernie, mediana wynosi 43 lata, gdzie dla osób z udarem 71 lat. Dodatkowo występują zmienne odstające w grupie ludzi z udarem. Sam rozkład zmiennej zbliżony jest do rozkładu jednostajnego z "grubym" oszacowaniem.
plotting(data = df, x='hypertension', type='class_hist')
0
Można zauważyć, że znaczna większość pacjentów nie ma nadciśnienia. Co więcej, mimo braku zbalansowanego zbioru danych, więcej pacjentów z udarem mózgu nie miało przy tym nadciśnienia.
plotting(data = df, x='heart_disease', type='class_hist')
0
Podobna sytuacja jak dla zmiennej hypertension. Znaczna większość pacjentow nie ma chorób serca oraz mimo braku zbalansowania, większość osób z udarem także nie ma chorób serca.
plotting(data = df, x='ever_married', type='class_hist')
0
Zmienna ever_married jest zmienną kategoryczną, którą będzie trzeba zakodować w wartości (0,1). Większa część pacjentów była kiedykolwiek żonata oraz to w tej grupie udar mózgu występuje w większej ilości.
plotting(data = df, x='work_type', type='class_hist')
0
Zmienna work_type określająca rodzaj wykonywanej pracy jest także zmienną kategoryczną, którą zakodujemy w wartości numeryczne. Najwięcej pacjentów jest z kategorii 'private', następnie równiemiernie ilościowo rozłożeni między 'children', 'self_employed', 'govt_job'. Znikoma ilość nigdy nie pracowała. Największa ilość udarów znajduje się wśród sektora prywatnego oraz samo-zatrunionego. Niewielka część w sektorze publicznym. Pojedyncze wartości występują też w grupie dzieci, ale z poprzednich wniosków najprawdopodobniej są to obserwacje odstające. Można także spróbować zkoszykować tą zmienną i rozdzielić na sektor obciążony i nieobciążony względem ryzyka udaru mózgu.
plotting(data = df, x='Residence_type', type='class_hist')
0
Zmienna określa obszar zamieszkania z rozdziałem na wiejski i miejski. Widać, że podział jest równomierny oraz w każdym z podziałów jest zbliżona proporcja osób po udarze do osób bez udaru mózgu, więc zmienna ta niesie ze sobą najpewniej niewielką ilość informacji dla potencjalnego modelu.
plotting(data = df, x='avg_glucose_level', type='hist_box')
0
plotting(data = df, x='avg_glucose_level', type='class_hist_box')
0
Możemy zaobserwować, że najwięcej jest pacjentów ze średnim poziomem glukozy w przedziale 77.24 - I kwantyl, 114.09 - II kwantyl, a więc z poziomem prawidłowym lub stanem przedcukrzycowym. Wniosek ten jednak jest uogólniony, gdyż zakres ten jest nieco inny w zależności od pory czy też sposobu wykonywanego badania, a także wieku badanej osoby. Osoby z wyższym średnim poziomem glukozy we krwi są outlierami dla grupy bez udaru mózgu, natomiast dla pozostałej części osób mieszczą się one w IV kwantylu. Więc usunięcie tych obserwacji, lub też zmiana ich wartości spowodowałaby utratę być może znaczącej informacji. Zamiast tego spróbujemy wykorzystać tą zależność tworząc nową zmienną kategoryczną informującą o średnim zakresie glukozy (poziom prawidłowy, stan przedcukrzycowy, cukrzyca) oraz imputując nowe wartości, jak też transformując całą zmienną.
plotting(data = df, x='bmi', type='hist_box')
0
plotting(data = df, x='bmi', type='class_hist_box')
0
Rozkład zmiennej bmi przypomina rozkład normalny z cięższym prawym ogonem. Występują też outliery dla bmi > 47.5. Większość pacjentów zawiera się w przedziale 23.5 - 33.1 czyli jest z wagą prawidłową lub nadwagą I stopnia. O ile sama zmienna może nie nieść zbyt dużej informacji z uwagi, że rozkład wartości jest zbliżony dla osób z udarem, jak i bez udaru to niemniej możliwe są pewne istotne interakcje między pozostałymi zmiennymi jak chociażby średnim poziomem glukozy we krwi. W kontekście danych odstających również można stworzyć zmienną kategoryczną zgodnie z kategoriami wskaźników BMI.
Dodatkowo zmienna BMI zawiera braki danych, którymi należy się zająć. Istnieje na to wiele metod m.in.
Z racji, iż braków danych jest niewiele (zaledwie 4%) to wypełnimi ich statystyką opisową jako działaniem prostym, ale i skutecznym. Jak wspomniano, rozkład zmiennej jest zbliżony do normalnego więc nie ma większego znaczenia wybór między medianą, a średnią ale zważywszy na lekki ogon prawostronny wybrana zostanie mediana. Dla obserwacji odstających posłużymy się metodą IQR.
plotting(data = df, x='smoking_status', type='class_hist_box')
0
Ciężko wyciągnać tu szersze wnioski zważywszy na dość liczną grupę 'Unknown', Najliczniejsza grupa pacjentów nigdy nie paliła papierosów zakładając, że w grupie 'Unknown' wartości rozkładałyby się w miare równomiernie. Widać rozróżnienei dla osób z udarem, że w większości są to osoby niepalące lub palące sporadycznie. Sprawdźmy wiek osób dla grupy Unknown, być może w naturalny sposób wydzieli nam się grupa osób niepełnoletnich.
plotting(data=df[(df.age <= 16)], x='smoking_status', type='class_hist')
0
Widzimy, że dla osób niepełnoletnich wartość smoking_status jest w większości równa 'Unknown'. Należy coś z tą grupą osób zrobić, jednakże nie można po prostu usunąć tych obserwacji, gdyż udar mózgu dotyczyć może także dzieci. Z jednej strony można zaklasyfikować te osoby do grupy 'Never smoked', jednak jest wysoce prawdopodobne, że w rzeczywistości część osób może 'popalać'. Stworzymy więc nową kategorię 'Underage' dla tych pacjentów.
W inżynierii cech będziemy tworzyli nowe zmienne z już istniejących, zajmowali się brakiem danych czy też odstającymi obserwacjami. Do tego ostatniego skorzystamy z metody rozstępu międzyćwiartylkowego IQR. Jest to różnica między trzecim a pierwszym kwartylem. Ponieważ pomiędzy tymi kwartylami znajduje się z definicji 50% wszystkich obserwacji (położonych centralnie w rozkładzie), dlatego im większa szerokość rozstępu ćwiartkowego, tym większe zróżnicowanie cechy. Dodatkowo na końcu dokonamy balansowania naszego zbioru z wykorzystaniem 4 technik: SMOTE, SMOTE+Tomek, SMOTE ENN.
#głęboka kopia danych
df_clear = df.copy()
#imputacja mediany w miejsca NA
df_clear['bmi_no_nan'] = df_clear.bmi.fillna(df_clear.bmi.median())
#zmienna numeryczna bez outlierów
df_clear['bmi_no_outliers'] = df_clear.bmi_no_nan
df_clear = iqr_outliers(df_clear, 'bmi_no_outliers')
#log-transformacja bez outlierów
df_clear['log_bmi_no_outliers'] = np.log(df_clear.bmi_no_nan)
df_clear = iqr_outliers(df_clear, 'log_bmi_no_outliers')
#nowa zmienna kategoryczna zgodnie ze wskaźnikami BMI
df_clear.loc[(df_clear.bmi_no_outliers < 18.5), 'bmi_categorical'] = 'underweight'
df_clear.loc[(df_clear.bmi_no_outliers >= 18.5) & (df_clear.bmi_no_outliers < 24.9), 'bmi_categorical'] = 'normal'
df_clear.loc[(df_clear.bmi_no_outliers >= 24.9) & (df_clear.bmi_no_outliers < 29.9), 'bmi_categorical'] = 'overweight'
df_clear.loc[(df_clear.bmi_no_outliers >= 29.9) & (df_clear.bmi_no_outliers < 34.9), 'bmi_categorical'] = 'obese'
df_clear.loc[(df_clear.bmi_no_outliers >= 34.9), 'bmi_categorical'] = 'extremely obese'
#kodowanie medianą przedziałów
df_clear.loc[df_clear.bmi_categorical == 'underweight', 'bmi_categorical_encoded'] = int(df_clear.loc[df_clear.bmi_categorical == 'underweight','bmi_no_outliers'].median())
df_clear.loc[df_clear.bmi_categorical == 'normal', 'bmi_categorical_encoded'] = int(df_clear.loc[df_clear.bmi_categorical == 'normal','bmi_no_outliers'].median())
df_clear.loc[df_clear.bmi_categorical == 'overweight', 'bmi_categorical_encoded'] = int(df_clear.loc[df_clear.bmi_categorical == 'overweight','bmi_no_outliers'].median())
df_clear.loc[df_clear.bmi_categorical == 'obese', 'bmi_categorical_encoded'] = int(df_clear.loc[df_clear.bmi_categorical == 'obese','bmi_no_outliers'].median())
df_clear.loc[df_clear.bmi_categorical == 'extremely obese', 'bmi_categorical_encoded'] = int(df_clear.loc[df_clear.bmi_categorical == 'extremely obese','bmi_no_outliers'].median())
df_clear.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 5110 entries, 0 to 5109 Data columns (total 17 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 id 5110 non-null int64 1 gender 5110 non-null object 2 age 5110 non-null float64 3 hypertension 5110 non-null int64 4 heart_disease 5110 non-null int64 5 ever_married 5110 non-null object 6 work_type 5110 non-null object 7 Residence_type 5110 non-null object 8 avg_glucose_level 5110 non-null float64 9 bmi 4909 non-null float64 10 smoking_status 5110 non-null object 11 stroke 5110 non-null int64 12 bmi_no_nan 5110 non-null float64 13 bmi_no_outliers 5110 non-null float64 14 log_bmi_no_outliers 5110 non-null float64 15 bmi_categorical 5110 non-null object 16 bmi_categorical_encoded 5110 non-null float64 dtypes: float64(7), int64(4), object(6) memory usage: 678.8+ KB
for col, types in zip(['bmi_no_nan', 'bmi_categorical', 'bmi_no_outliers', 'log_bmi_no_outliers'], ['class_hist_box','class_hist','class_hist_box','class_hist_box']):
plotting(data = df_clear, x=col, type=types)
Możemy zaobserwować, że wypełnienie braków danych spowodowało peak w punkcie mediany co jest naturalnym efektem zastosowanej metody.
Dodatkowo, najbardziej liczna jest grupa osób z nadwagą i to tam również najliczniej występują osoby z udarem mózgu. W przypadku usunięcia outlierów mozemy zaobserwować pozbycie się ogona prawostronnego, natomiast zastosowanie log-transformacji wydaje się nie mieć znacznego wpływu poza przeskalowaniem co będzie zbędne z racji iż w późniejszym etapie będziemy dokonywać standaryzacji zmiennych, dlatego odrzucimy tą zmienną.
Co więcej, możemy zauważyć, że przypisanie granicznych wartości kwantylowych nie zlikwidowało danych odstających w grupie osób z udarem. W tej sytuacji możemy usunąć obserwacje specjalnie dla tej grupy, ale nie są to pojedyncze przypadki a i tak jest ona mało liczna. W takim razie, z racji dalszego lekkiego zaszumienia zmiennej przez te obserwacje, pozostaniemy przy zmiennej kategorycznej koszykującej do odpowiedniej grupy. Jest to dość powszechne rozwiązanie dla BMI, gdyż to powszechny indeks wagi ciała.
df_clear = df_clear.loc[df_clear.gender != 'Other',:]
df_clear.loc[df_clear.gender == 'Female', 'is_female'] = 1
df_clear.loc[df_clear.gender == 'Male', 'is_female'] = 0
df_clear.is_female = df_clear.is_female.astype('int')
df_clear.info()
<class 'pandas.core.frame.DataFrame'> Int64Index: 5109 entries, 0 to 5109 Data columns (total 18 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 id 5109 non-null int64 1 gender 5109 non-null object 2 age 5109 non-null float64 3 hypertension 5109 non-null int64 4 heart_disease 5109 non-null int64 5 ever_married 5109 non-null object 6 work_type 5109 non-null object 7 Residence_type 5109 non-null object 8 avg_glucose_level 5109 non-null float64 9 bmi 4908 non-null float64 10 smoking_status 5109 non-null object 11 stroke 5109 non-null int64 12 bmi_no_nan 5109 non-null float64 13 bmi_no_outliers 5109 non-null float64 14 log_bmi_no_outliers 5109 non-null float64 15 bmi_categorical 5109 non-null object 16 bmi_categorical_encoded 5109 non-null float64 17 is_female 5109 non-null int32 dtypes: float64(7), int32(1), int64(4), object(6) memory usage: 738.4+ KB
#tablica z pacjentami z udarem
df_stroke = df_clear[df_clear.stroke == 1].copy()
#pozbycie sie outlierow
df_stroke = iqr_outliers(df=df_stroke, col='age')
#ponowna konkatenacja zbioru
df_clear = pd.concat([df_clear[df_clear.stroke==0], df_stroke])
#plot
plotting(data=df_clear, x='age', type='class_hist_box')
#podsumowanie
df_clear.info()
<class 'pandas.core.frame.DataFrame'> Int64Index: 5109 entries, 249 to 248 Data columns (total 18 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 id 5109 non-null int64 1 gender 5109 non-null object 2 age 5109 non-null float64 3 hypertension 5109 non-null int64 4 heart_disease 5109 non-null int64 5 ever_married 5109 non-null object 6 work_type 5109 non-null object 7 Residence_type 5109 non-null object 8 avg_glucose_level 5109 non-null float64 9 bmi 4908 non-null float64 10 smoking_status 5109 non-null object 11 stroke 5109 non-null int64 12 bmi_no_nan 5109 non-null float64 13 bmi_no_outliers 5109 non-null float64 14 log_bmi_no_outliers 5109 non-null float64 15 bmi_categorical 5109 non-null object 16 bmi_categorical_encoded 5109 non-null float64 17 is_female 5109 non-null int32 dtypes: float64(7), int32(1), int64(4), object(6) memory usage: 738.4+ KB
df_clear.loc[df_clear.ever_married == 'Yes', 'is_married'] = 1
df_clear.loc[df_clear.ever_married == 'No', 'is_married'] = 0
df_clear.is_married = df_clear.is_married.astype('int')
df_clear.info()
<class 'pandas.core.frame.DataFrame'> Int64Index: 5109 entries, 249 to 248 Data columns (total 19 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 id 5109 non-null int64 1 gender 5109 non-null object 2 age 5109 non-null float64 3 hypertension 5109 non-null int64 4 heart_disease 5109 non-null int64 5 ever_married 5109 non-null object 6 work_type 5109 non-null object 7 Residence_type 5109 non-null object 8 avg_glucose_level 5109 non-null float64 9 bmi 4908 non-null float64 10 smoking_status 5109 non-null object 11 stroke 5109 non-null int64 12 bmi_no_nan 5109 non-null float64 13 bmi_no_outliers 5109 non-null float64 14 log_bmi_no_outliers 5109 non-null float64 15 bmi_categorical 5109 non-null object 16 bmi_categorical_encoded 5109 non-null float64 17 is_female 5109 non-null int32 18 is_married 5109 non-null int32 dtypes: float64(7), int32(2), int64(4), object(6) memory usage: 758.4+ KB
le = preprocessing.LabelEncoder()
le.fit(df_clear.work_type)
df_clear['work_type_encoded'] = le.transform(df_clear.work_type)
df_clear.info()
<class 'pandas.core.frame.DataFrame'> Int64Index: 5109 entries, 249 to 248 Data columns (total 20 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 id 5109 non-null int64 1 gender 5109 non-null object 2 age 5109 non-null float64 3 hypertension 5109 non-null int64 4 heart_disease 5109 non-null int64 5 ever_married 5109 non-null object 6 work_type 5109 non-null object 7 Residence_type 5109 non-null object 8 avg_glucose_level 5109 non-null float64 9 bmi 4908 non-null float64 10 smoking_status 5109 non-null object 11 stroke 5109 non-null int64 12 bmi_no_nan 5109 non-null float64 13 bmi_no_outliers 5109 non-null float64 14 log_bmi_no_outliers 5109 non-null float64 15 bmi_categorical 5109 non-null object 16 bmi_categorical_encoded 5109 non-null float64 17 is_female 5109 non-null int32 18 is_married 5109 non-null int32 19 work_type_encoded 5109 non-null int32 dtypes: float64(7), int32(3), int64(4), object(6) memory usage: 778.3+ KB
df_clear.loc[df_clear.Residence_type == 'Rural', 'is_rural'] = 1
df_clear.loc[df_clear.Residence_type == 'Urban', 'is_rural'] = 0
df_clear.is_rural = df_clear.is_rural.astype('int')
df_clear.info()
<class 'pandas.core.frame.DataFrame'> Int64Index: 5109 entries, 249 to 248 Data columns (total 21 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 id 5109 non-null int64 1 gender 5109 non-null object 2 age 5109 non-null float64 3 hypertension 5109 non-null int64 4 heart_disease 5109 non-null int64 5 ever_married 5109 non-null object 6 work_type 5109 non-null object 7 Residence_type 5109 non-null object 8 avg_glucose_level 5109 non-null float64 9 bmi 4908 non-null float64 10 smoking_status 5109 non-null object 11 stroke 5109 non-null int64 12 bmi_no_nan 5109 non-null float64 13 bmi_no_outliers 5109 non-null float64 14 log_bmi_no_outliers 5109 non-null float64 15 bmi_categorical 5109 non-null object 16 bmi_categorical_encoded 5109 non-null float64 17 is_female 5109 non-null int32 18 is_married 5109 non-null int32 19 work_type_encoded 5109 non-null int32 20 is_rural 5109 non-null int32 dtypes: float64(7), int32(4), int64(4), object(6) memory usage: 798.3+ KB
#nowa zmienna kategoryczna zgodnie ze wskaźnikami BMI
df_clear.loc[(df_clear.avg_glucose_level < 140), 'glucose_categorical'] = 'normal'
df_clear.loc[(df_clear.avg_glucose_level >= 140) & (df_clear.avg_glucose_level <= 199), 'glucose_categorical'] = 'prediabetes'
df_clear.loc[(df_clear.avg_glucose_level > 199), 'glucose_categorical'] = 'diabetes'
#kodowanie
le_gl = preprocessing.LabelEncoder()
le_gl.fit(df_clear.glucose_categorical)
df_clear['glucose_categorical_encoded'] = le_gl.transform(df_clear.glucose_categorical)
df_clear.info()
#zmienna numeryczna bez outlierów
df_clear['glucose_no_outliers'] = df_clear.avg_glucose_level
df_clear = iqr_outliers(df_clear, 'glucose_no_outliers')
#log-transformacja bez outlierów
df_clear['log_glucose_no_outliers'] = np.log(df_clear.avg_glucose_level)
df_clear = iqr_outliers(df_clear, 'log_glucose_no_outliers')
<class 'pandas.core.frame.DataFrame'> Int64Index: 5109 entries, 249 to 248 Data columns (total 23 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 id 5109 non-null int64 1 gender 5109 non-null object 2 age 5109 non-null float64 3 hypertension 5109 non-null int64 4 heart_disease 5109 non-null int64 5 ever_married 5109 non-null object 6 work_type 5109 non-null object 7 Residence_type 5109 non-null object 8 avg_glucose_level 5109 non-null float64 9 bmi 4908 non-null float64 10 smoking_status 5109 non-null object 11 stroke 5109 non-null int64 12 bmi_no_nan 5109 non-null float64 13 bmi_no_outliers 5109 non-null float64 14 log_bmi_no_outliers 5109 non-null float64 15 bmi_categorical 5109 non-null object 16 bmi_categorical_encoded 5109 non-null float64 17 is_female 5109 non-null int32 18 is_married 5109 non-null int32 19 work_type_encoded 5109 non-null int32 20 is_rural 5109 non-null int32 21 glucose_categorical 5109 non-null object 22 glucose_categorical_encoded 5109 non-null int32 dtypes: float64(7), int32(5), int64(4), object(7) memory usage: 858.2+ KB
for col, types in zip(['glucose_categorical', 'glucose_no_outliers', 'log_glucose_no_outliers'], ['class_hist','class_hist_box','class_hist_box']):
plotting(data = df_clear, x=col, type=types)
Najbardziej liczna jest grupa osób z prawidłowym kategorią średniego poziomu glukozy we krwi i to tam również najliczniej występują osoby z udarem mózgu (zaliczono tam także I stopień otyłości).
W kontekście kategoryzacji zmiennej otrzymano raczej spadek informacji, gdyż została jedna najliczniejsza grupa osób z udarem oraz pozostałe z podobną licznością. Co więcej, jak wspomniano, ciężko jednoznacznie wpisać w kategorie gdyż nie znamy właściwości pomiaru cukru, jak też dodatkowych chorób pacjentów takich jak cukrzyca.
Widzimy także, że metoda IQR z imputacją zamiast usunięcia danych dalej zachowała pewne dane odstające. Jednakże przy skorzystaniu z log-transformacji, są to nieliczne wartości mocno zbliżone do wartości górnej. Mając to na uwadze, jak też fakt spadku informacji po kategoryzacji ograniczymy się do tej zmiennej.
df_clear['smoking_status_new'] = df_clear.smoking_status
df_clear.loc[df_clear.age <= 21, 'smoking_status_new'] = 'underage'
print("Grupa 'Unknown' stanowi %.2f procent obserwacji całościowego zbioru" % (len(df_clear[df_clear.smoking_status_new == 'underage'])/len(df_clear)))
plotting(data=df_clear, x='smoking_status_new', type='class_hist')
df_clear.info()
Grupa 'Unknown' stanowi 0.21 procent obserwacji całościowego zbioru
<class 'pandas.core.frame.DataFrame'> Int64Index: 5109 entries, 249 to 248 Data columns (total 26 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 id 5109 non-null int64 1 gender 5109 non-null object 2 age 5109 non-null float64 3 hypertension 5109 non-null int64 4 heart_disease 5109 non-null int64 5 ever_married 5109 non-null object 6 work_type 5109 non-null object 7 Residence_type 5109 non-null object 8 avg_glucose_level 5109 non-null float64 9 bmi 4908 non-null float64 10 smoking_status 5109 non-null object 11 stroke 5109 non-null int64 12 bmi_no_nan 5109 non-null float64 13 bmi_no_outliers 5109 non-null float64 14 log_bmi_no_outliers 5109 non-null float64 15 bmi_categorical 5109 non-null object 16 bmi_categorical_encoded 5109 non-null float64 17 is_female 5109 non-null int32 18 is_married 5109 non-null int32 19 work_type_encoded 5109 non-null int32 20 is_rural 5109 non-null int32 21 glucose_categorical 5109 non-null object 22 glucose_categorical_encoded 5109 non-null int32 23 glucose_no_outliers 5109 non-null float64 24 log_glucose_no_outliers 5109 non-null float64 25 smoking_status_new 5109 non-null object dtypes: float64(9), int32(5), int64(4), object(8) memory usage: 977.9+ KB
W dalszym ciągu grupa Unknown jest dość liczna, ale wydobyliśmy za to bardziej zróżnicowany podział i pozostawimy to w ten sposób. Dzięki temu możemy dostrzec, że w prawie każdej grupie mamy równy proporcję osób z udarem, jak też bez udaru. Nie jest to zbyt informacyjna zmienna wyszczególniająca informację na temat pacjentów z udarem. Moglibyśmy tutaj wyodrebnić grupę 'underage', a pozostałe zagregować otrzymując zmienną binarną, jednak będzie ona opierała się tak naprawdę tylko na wieku, a nie na rzeczywistym statusie palenia. Odrzucimy więc tą zmienną w dalszej analizie.
# exclude categorical data
plotting(data=df_clear.loc[:,~df_clear.columns.isin(
['id', 'bmi_categorical_encoded','work_type_encoded','glucose_categorical_encoded'])], x='none', type='heatmap')
0
Widoczna jest słaba korelacja między zmiennymi oraz między zmiennymi, a zmienną celu. Największe wartości widoczne są przy tych samych zmiennych, ale po przeprowadzeniu transformacji lub kodowania co jest rezultatem oczywistym - jedynie zmienna is_married koreluje ze zmienną age ale nie będziemy jednak wykluczać żadnej z nich. Nie występuje więc tutaj problem współliniowości. Ciężko także na tej podstawie dokonać selekcji zmiennych, więc raczej skorzystamy w późniejszym etapie z regularyzacji LASSO.
Bardzo często wejściowe zmienne są interaktywne ze sobą w sposób nieliniowy. Takie interkacje mogą zostać później wykryte i wykorzystane przez model przy etapie selekcji zmiennych. Co więcej, podniesienie niektórych zmiennych do pewnej potęgi może pomóc w wyodrębnieniu istotnych relacji między tą zmienną, a zmienną celu. Przeważnie wystarczy podniesienie do stopnia drugiego (ewentualnie trzeciego), gdyż wielomiany stopnia czwartego i wyższego są nadmiernie elastyczne, mają nieregularne kształty i mogą prowadzić do overfittingu. My wykorzystamy cechy wielomianowe stopnia drugiego.
#get PolynomialFeatures transfomer
poly = preprocessing.PolynomialFeatures(2)
#select only apriori chosen variables
X_poly = df_clear.select_dtypes(['int32','int64','float32','float64'])
X_poly = X_poly.loc[:,~X_poly.columns.isin(['id', 'bmi_no_nan','bmi','avg_glucose_level','glucose_no_outliers', 'log_bmi_no_outliers', 'bmi_no_outliers','glucose_categorical_encoded'])]
print("Ramka wejściowa")
X_poly.info()
print("\nRamka wyjściowa po PolyFeatures transformacji")
df_poly = pd.DataFrame(poly.fit_transform(X_poly), columns = poly.get_feature_names_out(X_poly.columns))
df_poly.info()
Ramka wejściowa <class 'pandas.core.frame.DataFrame'> Int64Index: 5109 entries, 249 to 248 Data columns (total 10 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 age 5109 non-null float64 1 hypertension 5109 non-null int64 2 heart_disease 5109 non-null int64 3 stroke 5109 non-null int64 4 bmi_categorical_encoded 5109 non-null float64 5 is_female 5109 non-null int32 6 is_married 5109 non-null int32 7 work_type_encoded 5109 non-null int32 8 is_rural 5109 non-null int32 9 log_glucose_no_outliers 5109 non-null float64 dtypes: float64(3), int32(4), int64(3) memory usage: 359.2 KB Ramka wyjściowa po PolyFeatures transformacji <class 'pandas.core.frame.DataFrame'> RangeIndex: 5109 entries, 0 to 5108 Data columns (total 66 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 1 5109 non-null float64 1 age 5109 non-null float64 2 hypertension 5109 non-null float64 3 heart_disease 5109 non-null float64 4 stroke 5109 non-null float64 5 bmi_categorical_encoded 5109 non-null float64 6 is_female 5109 non-null float64 7 is_married 5109 non-null float64 8 work_type_encoded 5109 non-null float64 9 is_rural 5109 non-null float64 10 log_glucose_no_outliers 5109 non-null float64 11 age^2 5109 non-null float64 12 age hypertension 5109 non-null float64 13 age heart_disease 5109 non-null float64 14 age stroke 5109 non-null float64 15 age bmi_categorical_encoded 5109 non-null float64 16 age is_female 5109 non-null float64 17 age is_married 5109 non-null float64 18 age work_type_encoded 5109 non-null float64 19 age is_rural 5109 non-null float64 20 age log_glucose_no_outliers 5109 non-null float64 21 hypertension^2 5109 non-null float64 22 hypertension heart_disease 5109 non-null float64 23 hypertension stroke 5109 non-null float64 24 hypertension bmi_categorical_encoded 5109 non-null float64 25 hypertension is_female 5109 non-null float64 26 hypertension is_married 5109 non-null float64 27 hypertension work_type_encoded 5109 non-null float64 28 hypertension is_rural 5109 non-null float64 29 hypertension log_glucose_no_outliers 5109 non-null float64 30 heart_disease^2 5109 non-null float64 31 heart_disease stroke 5109 non-null float64 32 heart_disease bmi_categorical_encoded 5109 non-null float64 33 heart_disease is_female 5109 non-null float64 34 heart_disease is_married 5109 non-null float64 35 heart_disease work_type_encoded 5109 non-null float64 36 heart_disease is_rural 5109 non-null float64 37 heart_disease log_glucose_no_outliers 5109 non-null float64 38 stroke^2 5109 non-null float64 39 stroke bmi_categorical_encoded 5109 non-null float64 40 stroke is_female 5109 non-null float64 41 stroke is_married 5109 non-null float64 42 stroke work_type_encoded 5109 non-null float64 43 stroke is_rural 5109 non-null float64 44 stroke log_glucose_no_outliers 5109 non-null float64 45 bmi_categorical_encoded^2 5109 non-null float64 46 bmi_categorical_encoded is_female 5109 non-null float64 47 bmi_categorical_encoded is_married 5109 non-null float64 48 bmi_categorical_encoded work_type_encoded 5109 non-null float64 49 bmi_categorical_encoded is_rural 5109 non-null float64 50 bmi_categorical_encoded log_glucose_no_outliers 5109 non-null float64 51 is_female^2 5109 non-null float64 52 is_female is_married 5109 non-null float64 53 is_female work_type_encoded 5109 non-null float64 54 is_female is_rural 5109 non-null float64 55 is_female log_glucose_no_outliers 5109 non-null float64 56 is_married^2 5109 non-null float64 57 is_married work_type_encoded 5109 non-null float64 58 is_married is_rural 5109 non-null float64 59 is_married log_glucose_no_outliers 5109 non-null float64 60 work_type_encoded^2 5109 non-null float64 61 work_type_encoded is_rural 5109 non-null float64 62 work_type_encoded log_glucose_no_outliers 5109 non-null float64 63 is_rural^2 5109 non-null float64 64 is_rural log_glucose_no_outliers 5109 non-null float64 65 log_glucose_no_outliers^2 5109 non-null float64 dtypes: float64(66) memory usage: 2.6 MB
Jak widzimy uzyskaliśmy wszystkie możliwe kombinacje zmiennych wielomianowych 2-giego stopnia. Wykluczymy jednak na start zmienne, których z pewnością nie wykorzystamy takie jak bias: '1' oraz zmienne binarne podniesione do kwadratu, jak też zmienne z interakcją 'stroke' w celu uniknięcia tzw. data leakage.
col_drop = list(df_poly.filter(regex='stroke'))[1:] + ['1','hypertension^2','heart_disease^2','is_female^2','is_married^2','is_rural^2']
df_poly = df_poly[df_poly.columns.drop(col_drop)]
cat_features = ['is_rural','is_married', 'is_female', 'hypertension', 'heart_disease', 'work_type_encoded']
col_filter = set(list(df_poly.filter(regex='age')) + list(df_poly.filter(regex='log_glucose_no_outliers')) + cat_features + ['bmi_categorical_encoded','stroke'])
df_poly = df_poly[col_filter]
df_poly.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 5109 entries, 0 to 5108 Data columns (total 27 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 heart_disease log_glucose_no_outliers 5109 non-null float64 1 is_rural 5109 non-null float64 2 work_type_encoded 5109 non-null float64 3 age hypertension 5109 non-null float64 4 is_married log_glucose_no_outliers 5109 non-null float64 5 is_female 5109 non-null float64 6 age bmi_categorical_encoded 5109 non-null float64 7 age is_female 5109 non-null float64 8 stroke 5109 non-null float64 9 age^2 5109 non-null float64 10 age log_glucose_no_outliers 5109 non-null float64 11 age 5109 non-null float64 12 is_female log_glucose_no_outliers 5109 non-null float64 13 is_rural log_glucose_no_outliers 5109 non-null float64 14 bmi_categorical_encoded 5109 non-null float64 15 bmi_categorical_encoded log_glucose_no_outliers 5109 non-null float64 16 age is_rural 5109 non-null float64 17 log_glucose_no_outliers 5109 non-null float64 18 work_type_encoded log_glucose_no_outliers 5109 non-null float64 19 log_glucose_no_outliers^2 5109 non-null float64 20 is_married 5109 non-null float64 21 hypertension 5109 non-null float64 22 heart_disease 5109 non-null float64 23 age heart_disease 5109 non-null float64 24 age is_married 5109 non-null float64 25 age work_type_encoded 5109 non-null float64 26 hypertension log_glucose_no_outliers 5109 non-null float64 dtypes: float64(27) memory usage: 1.1 MB
Ważne jest, aby balansować dane tylko na zbiorze treningowym, gdyż w ogólności zbiór testowy jest zbiorem, który ma odzwierciedlać realne dane - a te zbalansowane nie są.
#wyciagniecie macierzy zmiennych X oraz wektora zmiennej celu Y
X = df_poly.loc[:,~df_poly.columns.isin(['stroke'])]
Y = df_poly.stroke
#train-test split, aby balansować tylko na zbiorze treningowym
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.33, random_state=123)
#for visualization
X_train['index'] = X_train.reset_index().index
df_imbalanced = X_train.copy()
df_imbalanced['stroke'] = y_train
#plot ukazujacy problem niezbalansowanych danych - klasa pacjentow z udarem jest prawie niewidoczna
plotting(data = df_imbalanced,x='index', y='age', type='scatter')
0
SMOTE wykorzystuje podejście algorytmiczne KNN tj. dla każdej obserwacji $x$ wybiera $k$ najbliższych sąsiadów i wybiera losowo jednego $x'$ z nich. Następnie różnica między koordynatami $x$ oraz $x'$ jest obliczana i przemnożona losowo przez liczby z przedziału $(0,1)$. Tak stworzona z różnic obserwacja $x''$ jest dodana do zbioru. Geometrycznie odwzorowuje to stworzenie nowej obserwacji na podstawie przesunięcia obserwacji bieżącej do jednego z sąsiadów.
# SMOTE
df_smote, X_smote, y_smote = balance_data(X=X_train, y=y_train, cat_features = cat_features, algorithm = 'SMOTE', strategy=1, k_neighbors = 7)
plotting(data = df_smote, x='index', y='age', type='scatter')
Liczba wszystkich obserwacji: 6530 Liczba obserwacji z klasy mniejszościowej: 3265 Procent klasy mniejszościowej do całości po zbalansowaniu: 0.5
0
SMOTE + Tomek jest techniką hybrydową, która łączy technikę oversamplingu algorytmu SMOTE z techniką undersamplingu algorytmu Tomek. Po zastosowaniu SMOTE'a istnieje możliwość nakładania się klastrów klas co może skutkować w overfittingu. Wówczas użyty zostaje dodatkowo algorytm Tomek, który polega na sparowaniu ze sobą próbek z różnych klas, a następnie obserwacji z klasy większościowej zostają usunięte w celu (niekiedy z obydwu) wyraźniejszej separacji klas.
# SMOTE Tomek
df_smtom, X_smtom, y_smtom = balance_data(X=X_train, y=y_train, cat_features = cat_features, algorithm='SMOTETOMEK', strategy = 1, k_neighbors=7)
plotting(data = df_smtom, x='index', y='age', type='scatter')
Liczba wszystkich obserwacji: 6432 Liczba obserwacji z klasy mniejszościowej: 3216 Procent klasy mniejszościowej do całości po zbalansowaniu: 0.5
0
SMOTE + ENN jest techniką hybrydową, która polega na dodatkowym usunięciu pewnej liczby obserwacji z przestrzeni cech. Gdy zbiór zostanie poddany oversamplingu przez SMOTE to następnie zostaje zastosowana technika ENN, która jest techniką undersamplingu. ENN (Edited Nearest Neighbor) polega na odnalezieniu najbliższych K-sąsiadów każdej z obserwacji. Następnie sprawdza czy klasa większościowa tych sąsiadów jest zgodna z klasą danej obserwacji, jeśli nie to zarówno obserwacja jak i sąsiedzi zostają usunięci. Dzięki takiemu podejściu oddzielenie klas od siebie jest wyraźniejsze, aniżeli w pozostałych przypadkach.
# SMOTE ENN
df_smenn, X_smenn, y_smenn = balance_data(X=X_train, y=y_train, cat_features = cat_features, algorithm='SMOTEENN', strategy = 1)
plotting(data = df_smenn, x='index', y='age',type='scatter')
Liczba wszystkich obserwacji: 6320 Liczba obserwacji z klasy mniejszościowej: 3055 Procent klasy mniejszościowej do całości po zbalansowaniu: 0.4833860759493671
0
Widzimy, że każdy z algorytmów balansuje w specyficzny dla siebie sposób generuje sztuczne obserwacje dla klasy mniejszościowej. W algorytmie SMOTE ENN widać dodatkowo zmalenie liczby obserwacji za sprawą wyraźnej techniki undersamplingu. Przeprowadzimy trenowanie modeli na każdym ze zbiorze dla każdego z modeli, aby mieć szerszy widok na rezultaty. Co więcej dla każdego z nich przeprowadzimy także regularyzację LASSO.
W ogólności dla oceny rozwiązań problemów klasyfikacyjnych korzysta się z macierzy pomyłek.
Zestawia ona rzeczywiste wartości klas z wartościami klas zamodelowanymi. W przypadku binarnym można to przedstawić w postaci TP czyli przypadki zaklasyfikowane prawidłowo do klasy prawdziwej, TN czyli przypadki zaklasyfikowane prawidłowo do klasy fałszywej, FP czyli przypadki zaklasyfikowane nieprawidłowo do klasy prawdziwej oraz FN czyli przypadki zaklasyfikowane nieprawidłowo do klasy fałszywej. Wartości FP są tzw. błędem I typu, czyli zaklasyfikowaniem czegoś jako prawdziwe, gdy w rzeczywistości takie nie jest. Z kolei wartości FN są błędem II typu, czyli zaklasyfikowaniem czegoś jako nieprawdziwe, gdy w rzeczywistości prawdziwe jest.
Z takiego zestawienia najpopularniejsze są metryki:
W naszym przypadku chcemy się skupić na minimalizacji błędu II typu kosztem błędu typu I. Wynika to z faktu, iż w przypadku popełnienia błędu I typu czyli nieprawdziwej predykcji udaru mózgu mogą zostać wykonane niepotrzebne badania oraz poniesione środki finansowe, natomiast w przypadek błędu II typu niesie poważne konsekwencje zdrowotne a w najgorszym przypadku śmierć.
Nie wystarczy jednak minimalizacja FN, ponieważ w takim przypadku każdy kolejny pacjent byłby klasyfikowany jako osoba z udarem mózgu. Potrzebne jest pewne wyważenie między FP, a FN dlatego kierować się będziemy metrykami Precision - czyli ile osób zaklasyfikowanych z udarem rzeczywiście ma udar (wysoka wartość precyzji implikuje niską wartość FP) oraz Recall (wysoka wartość recall implikuje niską wartość FN), jak też ich średnią harmoniczną $F_1$.
Podsumowując, dla klasy pacjentów z udarem bardziej istotna jest metryka Recall, która informuje nas jaka jest proporcja osób zaklasyfikowanych z udarem do wszystkich osób, które zaklasyfikowane być powinny. Oczekujemy wysoką wartość, gdyż oznaczać to będzie że prawidłowo wyszukujemy osoby z udarem.
Dla pacjentów bez udaru bardziej patrzeć będziemy na metrykę Precision, która z kolei mówi nam ile zaklasyfikowanych jako osoby bez udaru rzeczywiście udaru nie mają. Wysoka wartość będzie oznaczać, że precyzyjnie wyszukujemy osoby bez udaru.
W naszym zestawie danych występuje wiele zmiennych po wprowadzeniu polynomial features. O ile w trakcie preprocessingu danych nie usunęliśmy ich zbyt wiele i wciąż liczba zmiennych jest o wiele mniejsza od liczby obserwacji to wciąż nadmierna ich ilość prowadzić może do zjawiska overfittingu, wpływać na wyjaśnialność, czy też wprowadzać niepotrzebny szum do modelu.
Dlatego też chcemy tą liczbę zredukować poprzez wybranie naijstotniejszych zmiennych (chociaż możliwe są inne sposoby jak redukcja wymiarowości przykładowo poprzez PCA).
Wytrenujemy 3 modele:
Regresje logistyczną
SVM czyli maszyne wektorów nośnych.
RandomForest
Każdy z modeli zostanie wytrenowany po wcześniejszej selekcji cech z wykorzystaniem regresji logistycznej z regularyzacją LASSO. Regularyzacja ta zeruje pewne współczynniki i zostanie ona wytrenowana przy wykorzystaniu kilku parametrów regularyzujących spośród których wybrany zostanie najlepszy dla metryki ROC AUC.
Po wybraniu tych zmiennych, wytrenujemy na tym wyselekcjowanym zbiorze wszystkie modele wraz z przeszukaniem siatki hiperparametrów każdego z nich maksymalizując wynik ROC AUC co pomoże nam w doborze ostatecznego modelu.
W ostatnim kroku przejdziemy do doboru thresholdu $p$, czyli wartości prawdopodobieństwa od którego klasyfikować będziemy pacjentów z udarem. Dla regresji logistycznej oraz RFC jest to naturalne, gdyż metody te są z natury probabilistyczne a ich algorytmy zwracają wartości prawdopodobieństw klasyfikacji. Dla SVM możemy ustawić wartość hiperparametru probability=True, dzięki czemu jesteśmy w stanie pozyskać owe prawdopodobieństwa.
Całość odbywa się za scenami przy wykorzystaniu regresji logistycznej na danych wynikach SVM (algorytm Platta). Technika ta nie jest do końca wiarygodna i ma zastrzeżenia, mianowicie:
W kontekście zastrzeżeń ustawiając własny threshold będziemy to w pewien sposób kontrolować tym samym ograniczając efektywność działania samego algorytmu SVM, a skupiając się na skuteczności działania techniki Platta w estymacji tych prawdopodobieństw.
report = pd.DataFrame()
lr_base = LogisticRegressionWithThreshold(penalty='l1', solver='saga',max_iter=10000, random_state = 123)
for X, y, name in zip([X_smote, X_smtom, X_smenn],[y_smote,y_smtom, y_smenn], ['SMOTE', 'SMOTE+TOMEK','SMOTE+ENN']):
print('\n'+name +'\n')
columns_to_standardize = list(set(list(X.filter(regex='age')) + list(X.filter(regex='log_glucose_no_outliers'))))
results = fit_pred_score(X, y, X_test, y_test, lr_base, dataset_name = name, model_name = 'Base LR', scaler = True, ohe=True, columns_to_standardize = columns_to_standardize, cat_features = cat_features, visualize_test = True, visualize_train = True)
report = pd.concat([report, results])
report.reset_index(inplace=True, drop=True)
report.sort_values(by=['model_name','metric','dataset_name'])
SMOTE
SMOTE+TOMEK
SMOTE+ENN
| No stroke | Stroke | accuracy | macro avg | weighted avg | metric | dataset_name | model_name | split | |
|---|---|---|---|---|---|---|---|---|---|
| 2 | 0.840624 | 0.840233 | 0.840429 | 0.840429 | 0.840429 | f1-score | SMOTE | Base LR | train |
| 6 | 0.908547 | 0.246575 | 0.836892 | 0.577561 | 0.872818 | f1-score | SMOTE | Base LR | test |
| 18 | 0.821209 | 0.820244 | 0.820728 | 0.820727 | 0.820743 | f1-score | SMOTE+ENN | Base LR | train |
| 22 | 0.886301 | 0.265487 | 0.803084 | 0.575894 | 0.852794 | f1-score | SMOTE+ENN | Base LR | test |
| 10 | 0.842220 | 0.842171 | 0.842195 | 0.842195 | 0.842195 | f1-score | SMOTE+TOMEK | Base LR | train |
| 14 | 0.908364 | 0.258760 | 0.836892 | 0.583562 | 0.873302 | f1-score | SMOTE+TOMEK | Base LR | test |
| 0 | 0.839597 | 0.841265 | 0.840429 | 0.840431 | 0.840431 | precision | SMOTE | Base LR | train |
| 4 | 0.967422 | 0.164234 | 0.836892 | 0.565828 | 0.924071 | precision | SMOTE | Base LR | test |
| 16 | 0.847005 | 0.795874 | 0.820728 | 0.821440 | 0.822289 | precision | SMOTE+ENN | Base LR | train |
| 20 | 0.976604 | 0.166205 | 0.803084 | 0.571404 | 0.932863 | precision | SMOTE+ENN | Base LR | test |
| 8 | 0.842089 | 0.842302 | 0.842195 | 0.842195 | 0.842195 | precision | SMOTE+TOMEK | Base LR | train |
| 12 | 0.969417 | 0.171429 | 0.836892 | 0.570423 | 0.926346 | precision | SMOTE+TOMEK | Base LR | test |
| 1 | 0.841654 | 0.839204 | 0.840429 | 0.840429 | 0.840429 | recall | SMOTE | Base LR | train |
| 5 | 0.856426 | 0.494505 | 0.836892 | 0.675466 | 0.836892 | recall | SMOTE | Base LR | test |
| 17 | 0.796937 | 0.846154 | 0.820728 | 0.821546 | 0.820728 | recall | SMOTE+ENN | Base LR | train |
| 21 | 0.811285 | 0.659341 | 0.803084 | 0.735313 | 0.803084 | recall | SMOTE+ENN | Base LR | test |
| 9 | 0.842351 | 0.842040 | 0.842195 | 0.842195 | 0.842195 | recall | SMOTE+TOMEK | Base LR | train |
| 13 | 0.854545 | 0.527473 | 0.836892 | 0.691009 | 0.836892 | recall | SMOTE+TOMEK | Base LR | test |
| 3 | 3265.000000 | 3265.000000 | 0.840429 | 6530.000000 | 6530.000000 | support | SMOTE | Base LR | train |
| 7 | 1595.000000 | 91.000000 | 0.836892 | 1686.000000 | 1686.000000 | support | SMOTE | Base LR | test |
| 19 | 3265.000000 | 3055.000000 | 0.820728 | 6320.000000 | 6320.000000 | support | SMOTE+ENN | Base LR | train |
| 23 | 1595.000000 | 91.000000 | 0.803084 | 1686.000000 | 1686.000000 | support | SMOTE+ENN | Base LR | test |
| 11 | 3216.000000 | 3216.000000 | 0.842195 | 6432.000000 | 6432.000000 | support | SMOTE+TOMEK | Base LR | train |
| 15 | 1595.000000 | 91.000000 | 0.836892 | 1686.000000 | 1686.000000 | support | SMOTE+TOMEK | Base LR | test |
plot_datasets_scores = transform_report_to_plot(report[(report.model_name=='Base LR') & (report.split=='test')], x='dataset_name', cl='metric')
plotting(data = plot_datasets_scores, x='', y='value', cl='metric', type='lineplot')
print("Wykres zależności wartości wybranych metryk od techniki balansowania")
Wykres zależności wartości wybranych metryk od techniki balansowania
report_diff = prepare_diff_dataset(report, by='dataset_name')
plot_diff_scores = transform_report_to_plot(report_diff[(report_diff.model_name=='Base LR') & (report_diff.train_test_diffaccuracy.isna()==False)], x='dataset_name', cl='metric', train_test=True)
plotting(data = plot_diff_scores, x='', y='value', cl='metric', type='lineplot')
print("Wykres zależności różnic między wartościami zbioru treningowego i testowego wybranych metryk od techniki balansowania")
Wykres zależności różnic między wartościami zbioru treningowego i testowego wybranych metryk od techniki balansowania
Dla ogólnego spojrzenia spoglądamy na metrykę 'macro-avg' jako, że nie uwzględnia proporcji klas (pod kątem liczby wystąpień), a chcemy aby klasa 'Stroke' była tak samo ważna w ocenie modelu.
Tak też dla:
W ogólności widać, że model sprawuje się całkiem dobrze bez przeprowadzenia preprocessingu oraz doboru hiperparametrów. W naszym przypadku wybór metody balansowania danych raczej nie jest zbyt istotny, ale dobrano 3 możliwości w celu pokazania różnych technik jak też pewnego wpływu na ewentualne wyniki. Tutaj różnice są naprawdę niewielkie, ale z racji na wyższe wartości, jak też niższe różnice względem zbioru treningowego i testowego wybierzemy metodą hybrydową SMOTE+ENN.
X = X_smenn.copy()
y = y_smenn.copy()
dataset_name = 'SMOTE+ENN'
columns_to_standardize = list(set(list(X.filter(regex='age')) + list(X.filter(regex='log_glucose_no_outliers')))) + ['bmi_categorical_encoded'] #bmi jest skategoryzowana, ale medianą w celu zachowania pewnej informacji w wartości dlatego ją standaryzujemy
X_train_preprocessed, X_test_preprocessed, standard_scaler = preprocess_data(X, X_test, scaler = True, ohe=True, columns_to_standardize = columns_to_standardize, cat_features = cat_features)
Do selekcji zmiennych wykorzystamy metodę LASSO w tym przypadku będzie to regresja logistyczna z funkcją kary w metryce L1. Przeszukamy siatkę parametrów regulujących C (im niższa wartość tym bardziej restrykcyjna funkcja kary). Będziemy maksymalizować wynik ROC AUC, czyli maksymalizować ilość prawdziwych pozytywnych predykcji przy jednoczesnej minimalizacji negatywnych predykcji. Następnie z najlepszego estymatora wybierzemy odpowiednie zmienne przy użyciu modułu SelectFromModel.
lr_grid = {'C': [1, 0.01, 0.001, 0.0001, 0.00001]}
lr = LogisticRegressionWithThreshold(penalty='l1', solver='saga',max_iter=10000, random_state = 123)
lr_search = GridSearchCV(estimator = lr, param_grid = lr_grid, scoring = 'f1',
cv = 3, n_jobs = -1)
lr_search.fit(X_train_preprocessed,y)
GridSearchCV(cv=3,
estimator=LogisticRegressionWithThreshold(max_iter=10000,
penalty='l1',
random_state=123,
solver='saga'),
n_jobs=-1, param_grid={'C': [1, 0.01, 0.001, 0.0001, 1e-05]},
scoring='f1')In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. GridSearchCV(cv=3,
estimator=LogisticRegressionWithThreshold(max_iter=10000,
penalty='l1',
random_state=123,
solver='saga'),
n_jobs=-1, param_grid={'C': [1, 0.01, 0.001, 0.0001, 1e-05]},
scoring='f1')LogisticRegressionWithThreshold(max_iter=10000, penalty='l1', random_state=123,
solver='saga')LogisticRegressionWithThreshold(max_iter=10000, penalty='l1', random_state=123,
solver='saga')lr = lr_search.best_estimator_
selection_model = SelectFromModel(lr).fit(X_train_preprocessed, y)
coefs = {}
for coef, feature in zip(selection_model.estimator.coef_[0], selection_model.estimator.feature_names_in_):
coefs[feature] = coef
X_selected = selection_model.transform(X_train_preprocessed)
print("Wymiar początkowego zbioru danych: ", X_train_preprocessed.shape)
print("Wymiar zbioru po selekcji zmiennych: ", X_selected.shape)
coefs
Wymiar początkowego zbioru danych: (6320, 26) Wymiar zbioru po selekcji zmiennych: (6320, 26)
{'heart_disease log_glucose_no_outliers': 0.8344674630857023,
'is_rural': -0.3836178773412988,
'work_type_encoded': 0.06631646886871749,
'age hypertension': -0.5360502118541486,
'is_married log_glucose_no_outliers': -1.3696022295920762,
'is_female': -0.5486473753957081,
'age bmi_categorical_encoded': -1.894616593787902,
'age is_female': -0.28609178877921704,
'age^2': -1.7626312272903806,
'age log_glucose_no_outliers': 3.542685050407526,
'age': 3.29947886926461,
'is_female log_glucose_no_outliers': 0.5107077922181565,
'is_rural log_glucose_no_outliers': 0.3124455072056057,
'bmi_categorical_encoded': -0.27217306234852345,
'bmi_categorical_encoded log_glucose_no_outliers': 1.213534468528651,
'age is_rural': -0.24886293590718478,
'log_glucose_no_outliers': -3.2697571834283563,
'work_type_encoded log_glucose_no_outliers': 0.29329202605187144,
'log_glucose_no_outliers^2': 2.320938764383206,
'is_married': 4.113051514840922,
'hypertension': -1.8745757969323773,
'heart_disease': -1.5783786106522788,
'age heart_disease': -0.392978351567816,
'age is_married': -0.20652025298948645,
'age work_type_encoded': -0.3296512727776314,
'hypertension log_glucose_no_outliers': 1.044660005331124}
Widzimy, że żadna zmienna nie została "usuniętą" przez LASSO poprzez przypisanie zerowej wagi.
Ostatnim etapem będzie dostrojenie hiperparametrów dla każdego z modeli. W regresji logistycznej dobierzemy wartości parametru $C$ odpowiadającego za regularyzację. Dla SVM wypróbujemy inne wartości gamma, $C$ oraz jądro liniowe. Dla RFC stworzymy przestrzeń hiperparametrów uwzględniającą m.in. głębokość drzewa. Wykorzystamy do tego GridSearchCV z pakietu sklearn. Metryką, którą będziemy optymalizować będzie metryka ROC AUC, która odzwierciedla wartości przedstawione na powyższym wykresie. Im wyższa wartość AUC tym wyższy odsetek klsayfikacji prawdziwie pozytywnych przy niskim odsetku klasyfikacji fałszywie pozytywnych. Ta metryka pozwala dobrze zestawiać ze sobą modele oraz podobnie jak poprzednie obrane - poprawnie odzwierciedla nasz cel.
svm_grid = {'C': [1, 0.01, 0.0001],
'gamma': [1, 0.01, 0.0001],
'kernel': ['rbf', 'linear']}
rfc_grid = {
'bootstrap': [True],
'max_depth': [80, 100],
'max_features': [2, 3],
'min_samples_leaf': [3, 5],
'min_samples_split': [8,12],
'n_estimators': [100, 500, 1000]
}
# Create a based model
svm = SupportVectorMachineWithThreshold()
rfc = RandomForestClassifierWithThreshold(random_state = 123)
# Instantiate the grid search model
svm_search = GridSearchCV(estimator = svm, param_grid = svm_grid, scoring = 'f1',
cv = 3, n_jobs = -1, verbose = 2)
rfc_search = GridSearchCV(estimator = rfc, param_grid = rfc_grid, scoring = 'f1',
cv = 3, n_jobs = -1, verbose = 2)
for grid in [svm_search, rfc_search]:
grid.fit(X_selected,y)
print(grid.best_estimator_)
Fitting 3 folds for each of 18 candidates, totalling 54 fits
SupportVectorMachineWithThreshold(C=1, gamma=1)
Fitting 3 folds for each of 48 candidates, totalling 144 fits
RandomForestClassifierWithThreshold(max_depth=80, max_features=3,
min_samples_leaf=3, min_samples_split=8,
random_state=123)
Otrzymaliśmy więc 3 modele z najlepszymi parametrami pod względem metryki ROC AUC z przeszukiwanych siatek parametrów. Teraz zestawimy je ze sobą.
plt.figure(0).clf()
svm = SupportVectorMachineWithThreshold(gamma=1, C=1, kernel='rbf', probability=True, random_state = 123)
rfc = RandomForestClassifierWithThreshold(max_depth=80, max_features=3,
min_samples_leaf=3, min_samples_split=8,
n_estimators=500, random_state=123)
#fit logistic regression model and plot ROC curve
lr.fit(X_selected, y)
svm.fit(X_selected, y)
rfc.fit(X_selected, y)
RandomForestClassifierWithThreshold(max_depth=80, max_features=3,
min_samples_leaf=3, min_samples_split=8,
n_estimators=500, random_state=123)In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. RandomForestClassifierWithThreshold(max_depth=80, max_features=3,
min_samples_leaf=3, min_samples_split=8,
n_estimators=500, random_state=123)<Figure size 576x396 with 0 Axes>
y_pred_lr = lr.predict(X_test_preprocessed)
y_pred_svm = svm.predict(X_test_preprocessed)
y_pred_rfc = rfc.predict(X_test_preprocessed)
fpr_lr, tpr_lr, _ = metrics.roc_curve(y_test, y_pred_lr)
auc = round(metrics.roc_auc_score(y_test, y_pred_lr), 4)
plt.plot(fpr_lr,tpr_lr,label="Logistic Regression, AUC="+str(auc))
fpr_svm, tpr_svm, _ = metrics.roc_curve(y_test, y_pred_svm)
auc = round(metrics.roc_auc_score(y_test, y_pred_svm), 4)
plt.plot(fpr_svm,tpr_svm,label="Support Vector Machine, AUC="+str(auc))
fpr_rfc, tpr_rfc, _ = metrics.roc_curve(y_test, y_pred_rfc)
auc = round(metrics.roc_auc_score(y_test, y_pred_rfc), 4)
plt.plot(fpr_rfc,tpr_rfc,label="Random Forest Classifier, AUC="+str(auc))
#add legend
plt.legend()
plt.show()
Widzimy, że spośród najlepszych estymatorów dla wybranej przestrzeni parametrów najlepszy wynik osiąga regresja logistyczna. Z kolei SVM oraz RFC sprawują się niewiele lepiej jak bazowy estymator czyli predykujący z prawdopodobieństwem 50%. To tak jakbyśmy typowali czy ktoś ma udar rzutem monetą - SVM i RFC sprawują się tylko niewiele lepiej dla tego scenariusza. Ograniczymy się więc do regresji logistycznej i sprawdzimy czy wybór thresholda $p$ wpłynie na nasz wynik.
for p in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]:
print('\nWartość p='+str(p)+' dla modelu Regresji Logistycznej\n')
model_name = 'Logistic Regression with p='+str(p)
results = fit_pred_score(X_selected, y, X_test_preprocessed, y_test, lr, dataset_name = dataset_name, model_name = model_name, scaler = False, ohe=False,
visualize_test = True, visualize_train = True, threshold = p)
report = pd.concat([report, results])
report.reset_index(inplace=True, drop=True)
report.sort_values(by=['model_name','metric','dataset_name'])
Wartość p=0.1 dla modelu Regresji Logistycznej
Wartość p=0.2 dla modelu Regresji Logistycznej
Wartość p=0.3 dla modelu Regresji Logistycznej
Wartość p=0.4 dla modelu Regresji Logistycznej
Wartość p=0.5 dla modelu Regresji Logistycznej
Wartość p=0.6 dla modelu Regresji Logistycznej
Wartość p=0.7 dla modelu Regresji Logistycznej
Wartość p=0.8 dla modelu Regresji Logistycznej
Wartość p=0.9 dla modelu Regresji Logistycznej
| No stroke | Stroke | accuracy | macro avg | weighted avg | metric | dataset_name | model_name | split | |
|---|---|---|---|---|---|---|---|---|---|
| 2 | 0.840624 | 0.840233 | 0.840429 | 0.840429 | 0.840429 | f1-score | SMOTE | Base LR | train |
| 6 | 0.908547 | 0.246575 | 0.836892 | 0.577561 | 0.872818 | f1-score | SMOTE | Base LR | test |
| 18 | 0.821209 | 0.820244 | 0.820728 | 0.820727 | 0.820743 | f1-score | SMOTE+ENN | Base LR | train |
| 22 | 0.886301 | 0.265487 | 0.803084 | 0.575894 | 0.852794 | f1-score | SMOTE+ENN | Base LR | test |
| 10 | 0.842220 | 0.842171 | 0.842195 | 0.842195 | 0.842195 | f1-score | SMOTE+TOMEK | Base LR | train |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 92 | 0.948214 | 0.666667 | 0.947212 | 0.807440 | 0.933018 | precision | SMOTE+ENN | Logistic Regression with p=0.9 | test |
| 89 | 0.996937 | 0.284124 | 0.652373 | 0.640531 | 0.652373 | recall | SMOTE+ENN | Logistic Regression with p=0.9 | train |
| 93 | 0.998746 | 0.043956 | 0.947212 | 0.521351 | 0.947212 | recall | SMOTE+ENN | Logistic Regression with p=0.9 | test |
| 91 | 3265.000000 | 3055.000000 | 0.652373 | 6320.000000 | 6320.000000 | support | SMOTE+ENN | Logistic Regression with p=0.9 | train |
| 95 | 1595.000000 | 91.000000 | 0.947212 | 1686.000000 | 1686.000000 | support | SMOTE+ENN | Logistic Regression with p=0.9 | test |
96 rows × 9 columns
plot_datasets_scores = transform_report_to_plot(report[(report.model_name.str.contains('Logistic Regression with p')) & (report.split=='test')], x='model_name')
plotting(data = plot_datasets_scores, x='', y='value', cl='metric', type='lineplot')
0
report_diff_p = prepare_diff_dataset(report, by='model_name')
plot_diff_scores_p = transform_report_to_plot(report_diff_p[(report_diff_p.model_name.str.contains('Logistic Regression with p')) & (report_diff_p.train_test_diffaccuracy.isna()==False)], x='model_name', cl='metric', train_test=True)
plotting(data = plot_diff_scores_p, x='', y='value', cl='metric', type='lineplot')
print("Wykres zależności różnic między wartościami zbioru treningowego i testowego wybranych metryk od wartości p")
#plot_diff_scores_p
Wykres zależności różnic między wartościami zbioru treningowego i testowego wybranych metryk od wartości p
Wnioski dla:
wartości Precision:
wartości Recall:
Jak wspomnieliśmy na początku, chcemy maksymalizować wartość Precyzji dla klasy bez udaru (0) oraz Recall dla klasy z udarem (1). Dodatkowo wysokie wartości Precyzji dla klasy (1) oraz Recall dla klasy (0) będą atutem.
Od wartości $p=0.2$ dla LR możemy zaobserwować mocny spadek Recall dla klasy z udarem. W punkcie spadku dla klasy (0) wartość precyzji jest relatywnie wysoka, a dla klasy (1) następuje niewielki spadek względem Recall. Dodatkowo niskie wartości $p$ mają niskie różnice metryk między zbiorem treningowym a testowym. Są też zgodne z "wyczuciem dziedzinowym" tj. chcemy być bardziej uważni na możliwość wystąpienia udaru i raczej wystarczy nam już niewielkie prawdopodobieństwo wystąpienia, aby zaklasyfikować pacjenta. Co prawda wzrasta ryzyko nieprawidziwej klasyfikacji, ale w najgorszym przypadku zostaną poniesione niepotrzebne koszty (z takiej racji nie może być to też wartość zbyt niska, gdyż prawie każdego od razu będziemy klasyfikować jako osobę z udarem). Z drugiej strony zbyt wysoka wartość wymaganego prawdopodobieństwa może narazić zdrowie, a nawet życie ludzkie (widać to chociażby po drastycznym spadku Recall dla klasy 1 co oznacza, że mocno wzrosła liczba osób zaklasyfikowanych jako bez udaru, gdzie w rzeczywistości on wystąpił).
Interesują nas wartości $p \in (0.2, 0.3)$, ponieważ mają wysokie wartości metryk Precision dla klasy (0) oraz Recall dla klasy (1), które ustaliliśmy na początku. Jednocześnie są to momenty, w których wartość Precision dla (1) oraz Recall dla (0) zaczyna wzrastać wraz z $F_1$ score.
Wartość $p=0.1$ klasyfikowałaby prawie każdego jako pacjenta z udarem co jest niepożądanym efektem i do takich wniosków nie jest potrzebny żaden model. Z kolei w wartościach $p > 0.3$ widzimy mocny spadek Recall dla klasy 1, a więc tracimy pewność że wszyscy pacjencji z udarem, którzy winni być zaklasyfikowani rzeczywiście zostali przypisani do tej klasy.
Zostaje decyzja między wyższymi wartościami ustawionych apriori metryk, a pewną intuicją dziedzinową, gdyż dla $p=0.2$ również istnieje ryzyko klasyfikacji zbyt wielu "zdrowych" pacjentów. Sprawdźmy więc wartości AUC dla tych dwóch parametrów, które poinformują nas jaki scenariusz maksymalizuje TPR przy minimalnej ilości FPR.
plt.figure(0).clf()
y_pred_lr1 = lr.predict(X_test_preprocessed, threshold=0.1)
y_pred_lr2 = lr.predict(X_test_preprocessed, threshold=0.2)
y_pred_lr3 = lr.predict(X_test_preprocessed, threshold=0.3)
y_pred_lr4 = lr.predict(X_test_preprocessed, threshold=0.4)
fpr_lr1, tpr_lr1, _ = metrics.roc_curve(y_test, y_pred_lr1)
auc = round(metrics.roc_auc_score(y_test, y_pred_lr1), 4)
plt.plot(fpr_lr1,tpr_lr1,label="Logistic Regression with p=0.1, AUC="+str(auc))
fpr_lr2, tpr_lr2, _ = metrics.roc_curve(y_test, y_pred_lr2)
auc = round(metrics.roc_auc_score(y_test, y_pred_lr2), 4)
plt.plot(fpr_lr2,tpr_lr2,label="Logistic Regression with p=0.2, AUC="+str(auc))
fpr_lr3, tpr_lr3, _ = metrics.roc_curve(y_test, y_pred_lr3)
auc = round(metrics.roc_auc_score(y_test, y_pred_lr3), 4)
plt.plot(fpr_lr3,tpr_lr3,label="Logistic Regression with p=0.3, AUC="+str(auc))
fpr_lr4, tpr_lr4, _ = metrics.roc_curve(y_test, y_pred_lr4)
auc = round(metrics.roc_auc_score(y_test, y_pred_lr4), 4)
plt.plot(fpr_lr4,tpr_lr4,label="Logistic Regression with p=0.4, AUC="+str(auc))
fpr_lr, tpr_lr, _ = metrics.roc_curve(y_test, y_pred_lr)
auc = round(metrics.roc_auc_score(y_test, y_pred_lr), 4)
plt.plot(fpr_lr,tpr_lr,label="Logistic Regression with p=0.5, AUC="+str(auc))
#add legend
plt.legend()
plt.show()
Jak możemy zauważyć, potwierdziły się intuicje wypisane w poprzednim komentarzu.
Osiągnęliśmy wartość $AUC=0.7652$ dla $p=0.2$ co jest lepszym wynikiem od domyślnej wartości parametru $p=0.5$, dla której $AUC=0.7314$, więc jest pewna poprawa modelu. Widzimy również fakt, że wysokie wartości metryk wybranych apriori słusznie maksymalizują wartość AUC przy czym widocznie uwzględniają ryzyko klasyfikowania większości pacjentów jako osoby z udarem. Ewidentnie jest to dobry kierunek w przypadku niezbalansowanego zbioru jak nasz.
Dla $p=0.3$ mamy $AUC=0.7521$ co także nie jest złym wynikiem i po konsultacji ze specjalistą być może zostałby wybrany, gdyż przypisuje klasę z wyższym poziomem prawdopodobieństwa.
Widać również jak ze wzrostem parametru $p$ maleje wartość $AUC$.
Ostatecznie dobieramy wartość $p=0.2$ jako wartość naszego głównego modelu. Ostatnim etapie przygotujemy cały gotowy pipeline do predykcji udaru.
#Save final model
filename = 'final_lr_no_smoking_status_model.sav'
joblib.dump(lr, filename)
['final_lr_no_smoking_status_model.sav']
def preprocess_final(X_new: pd.DataFrame, standard_scaler):
"""
Function for preprocessing raw data.
It performs even unwanted operations because of the lack of time to prepare cleaning the code.
"""
df_clear = X_new
### BMI
df_clear['bmi_no_nan'] = 28.1
df_clear['bmi_no_outliers'] = df_clear.bmi_no_nan
df_clear = iqr_outliers(df_clear, 'bmi_no_outliers')
df_clear.loc[(df_clear.bmi_no_outliers < 18.5), 'bmi_categorical_encoded'] = 17.0
df_clear.loc[(df_clear.bmi_no_outliers >= 18.5) & (df_clear.bmi_no_outliers < 24.9), 'bmi_categorical_encoded'] = 22.0
df_clear.loc[(df_clear.bmi_no_outliers >= 24.9) & (df_clear.bmi_no_outliers < 29.9), 'bmi_categorical_encoded'] = 27.0
df_clear.loc[(df_clear.bmi_no_outliers >= 29.9) & (df_clear.bmi_no_outliers < 34.9), 'bmi_categorical_encoded'] = 32.0
df_clear.loc[(df_clear.bmi_no_outliers >= 34.9), 'bmi_categorical_encoded'] = 38.0
### GENDER
df_clear = df_clear.loc[df_clear.gender != 'Other',:]
df_clear.loc[df_clear.gender == 'Female', 'is_female'] = 1
df_clear.loc[df_clear.gender == 'Male', 'is_female'] = 0
df_clear.is_female = df_clear.is_female.astype('int')
### AGE
df_clear = iqr_outliers(df=df_clear, col='age')
### EVER_MARRIED
df_clear.loc[df_clear.ever_married == 'Yes', 'is_married'] = 1
df_clear.loc[df_clear.ever_married == 'No', 'is_married'] = 0
df_clear.is_married = df_clear.is_married.astype('int')
### WORK_TYPE
le = preprocessing.LabelEncoder()
le.fit(df_clear.work_type)
df_clear['work_type_encoded'] = le.transform(df_clear.work_type)
### RESIDENCE_TYPE
df_clear.loc[df_clear.Residence_type == 'Rural', 'is_rural'] = 1
df_clear.loc[df_clear.Residence_type == 'Urban', 'is_rural'] = 0
df_clear.is_rural = df_clear.is_rural.astype('int')
### AVG_GLUCOSE_LEVEL
df_clear['log_glucose_no_outliers'] = np.log(df_clear.avg_glucose_level)
df_clear = iqr_outliers(df_clear, 'log_glucose_no_outliers')
### POLYNOMIAL FEATURES
poly = preprocessing.PolynomialFeatures(2)
#select only apriori chosen variables
X_poly = df_clear.select_dtypes(['int32','int64','float32','float64'])
X_poly = X_poly.loc[:,~X_poly.columns.isin(['id','smoking_status', 'bmi_no_nan','bmi','avg_glucose_level','glucose_no_outliers', 'log_bmi_no_outliers', 'bmi_no_outliers','glucose_categorical_encoded'])]
df_poly = pd.DataFrame(poly.fit_transform(X_poly), columns = poly.get_feature_names_out(X_poly.columns))
df_poly = df_poly[['age','age bmi_categorical_encoded','age heart_disease','age hypertension','age is_female',
'age is_married','age is_rural','age log_glucose_no_outliers','age work_type_encoded',
'age^2','bmi_categorical_encoded','bmi_categorical_encoded log_glucose_no_outliers','heart_disease','heart_disease log_glucose_no_outliers',
'hypertension','hypertension log_glucose_no_outliers','is_female','is_female log_glucose_no_outliers','is_married','is_married log_glucose_no_outliers',
'is_rural','is_rural log_glucose_no_outliers','log_glucose_no_outliers','log_glucose_no_outliers^2','work_type_encoded','work_type_encoded log_glucose_no_outliers']]
cat_features = ['is_rural','is_married', 'is_female', 'hypertension', 'heart_disease', 'work_type_encoded']
columns_to_standardize = list(set(list(df_poly.filter(regex='age')) + list(df_poly.filter(regex='log_glucose_no_outliers')))) + ['bmi_categorical_encoded'] #bmi jest skategoryzowana, ale medianą w celu zachowania pewnej informacji w wartości dlatego ją standaryzujemy
### STANDARDIZE + OHE
X_new_preprocessed = preprocess_data(None, df_poly, scaler = True, ohe=True, columns_to_standardize = columns_to_standardize, cat_features = cat_features, test_only=True, sc=standard_scaler)
return X_new_preprocessed
#X_new
def predict_stroke(X : pd.DataFrame(), lr_model, standard_scaler):
"""
Function predicts the strokes based on the given features using trained model.
"""
#create model with the best parameters
model = lr_model
#preprocess data
X_preprocessed = preprocess_final(X, standard_scaler)
predicted_strokes = model.predict(X_preprocessed, threshold=0.3)
return predicted_strokes
#Load previously trained model
loaded_model = joblib.load(filename)
loaded_model
LogisticRegressionWithThreshold(C=1, max_iter=10000, penalty='l1',
random_state=123, solver='saga')In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. LogisticRegressionWithThreshold(C=1, max_iter=10000, penalty='l1',
random_state=123, solver='saga')#Predict strokes
predict_stroke(df.loc[0:5,:], loaded_model, standard_scaler)
array([0, 0, 0, 0, 0, 0])